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Fluid motion in a fuel tank produced during thrust oscillations can circulate sub-cooled 
hydrogen near the liquid-vapor interface resulting in increased condensation and ullage 
pressure collapse. The first objective of this study is to validate the capabilities of a 
Computational Fluid Dynamics (CFD) tool, CFD-ACE+, in modeling the fundamental 
interface transition physics occurring at the propellant surface. The second objective is to use 
the tool to assess the effects of thrust oscillations on surface dynamics. Our technical 
approach is to first verify the CFD code against known theoretical solutions, and then validate 
against existing experiments for small scale tanks and a range of transition regimes. A 2D 
axisymmetric, multi-phase model of gases, liquids, and solids is used to verify that CFD-ACE+ 
is capable of modeling fluid-structure interaction and system resonance in a typical thrust 
oscillation environment. Then, the 3D mode is studied with an assumed oscillatory body 
force to simulate the thrust oscillating effect. The study showed that CFD modeling can 
capture all of the transition physics from solid body motion to standing surface wave and to 
droplet ejection from liquid-gas interface. Unlike the analytical solutions established during 
the 1960’s, CFD modeling is not limited to the small amplitude regime. It can extend solutions 
to the nonlinear regime to determine the amplitude of surface waves after the onset of 
instability. The present simulation also demonstrated consistent trends from numerical 
experiments through variation of physical properties from low viscous fluid to high viscous 
fluids, and through variation of geometry and input forcing functions. A comparison of 
surface wave patterns under various forcing frequencies and amplitudes showed good 
agreement with experimental observations. It is concluded that thrust oscillations can cause 
droplet formation at the interface, which results in increased surface area and enhanced heat 
transfer between the liquid and gas phases as the ejected droplets travel well into the warmer 
gas region. 


Nomenclature 

LH2 = Liquid Hydrogen 

LOX = Liquid Oxygen 


I. Introduction 

T hrust oscillation of a solid rocket motor, also called resonant burning, is a phenomenon characterized by 
increased acceleration pulses that may be felt during the latter stages of first-stage powered flight. These 
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longitudinal forces may increase the loads experienced by the Ares I vehicle during flight, and may exceed 
allowable loads on various portions of the vehicle and on the crew. Fluid motion induced during thrust oscillation 
can circulate sub-cooled hydrogen near the liquid vapor interface resulting in increased condensation of ullage gas 
and corresponding pressure collapse. Conversely, liquid contact on hot tank walls caused by sloshing can result in 
rapid vaporization and subsequent rapid pressure rise. Understanding the dynamic and thermodynamic behavior 
and the response of cryogenic propellants of Liquid Hydrogen (LH2) and Liquid Oxygen (LOX) tanks under 
spacecraft flight operations (such as engine shut-down and restart, oscillatory thrust and oscillatory side forces) is 
extremely important during the vehicle design process. In addition, the proper orientation of propellant affects the 
performance of several systems including the tank venting system to ensure sufficient liquid engine start-up, feed 
system thermal conditioning, and stability conditions in feed lines prior to 2 nd stage ignition. The orientation control 
of the liquid surface is also crucial to minimize vehicle disturbances and avoid excessive demands on the attitude 
control system [1]. 

The liquid surface instability of a propellant tank due to thrust oscillation can be investigated through several 
different approaches: 

Analytical Solutions: A representative example of an analytical solution is the well-known Mathieu’s 

Equation [2] for a cylindrical tank filled with liquid and subjected to a vertical oscillatory driving force. The 
solution provides important insight into the physics and dependencies on various variables. However, these 
solutions are only obtainable under restricted conditions and for simple geometries. Mathieu’s solution also applies 
only to the identification of instability conditions without providing any further insight into the dynamics after the 
onset of the instability. 

Experimental Testing: Experiments provide the most reliable and realistic information. However, experiments 
can be expensive and time consuming. This is especially true for full scale cryogenic testing. Normally, only small 
scale model and room temperature data are available. It is difficult to obtain measurements for realistic flight 
conditions. Only limited data (if any) exist for cryogenic fluids under thrust oscillation during flight. 

Computational Modeling: The computational modeling approach of discretizing and numerically solving a 

relevant set of fluid and solid physics laws (the multi-physics CFD approach) has no geometrical limitation and can 
be applied to realistic full scale propellant tank geometries under flight conditions. It is not limited to the small 
amplitude regime, and can extend beyond the analytical solutions to the nonlinear regime to determine the amplitude 
of surface waves after the onset of instability. This is significant in that the phenomena such as droplet formation 
and ejection from the surface, which are due to nonlinear instability, are important drivers leading to potential ullage 
collapse. Computational modeling can also be applied to a range of fluids such as water, LH2, LOX, or Helium. It 
can bridge the analytical solutions and experiments and become a powerful design tool after systematic validation. 

The purpose of this study is to investigate if a Computational Fluid Dynamics tool is capable of accurately 
modeling the typical transition phenomena acting on the propellant surface observed in experiments: transitioning 
from solid body motion to standing surface waves to droplets under increased forcing. We also assess how CFD 
can become a reliable tool in assisting in the design and planning of an integrated sloshing test proposed at NASA 
MSFC facilities. Our technical approach is to first verify a CFD program, CFD-ACE+, against known theoretical 
solutions, and then validate existing experiments for small scale tanks across a range of transition regimes. The 
verified and validated code is then applied to full scale propellant tanks to study the wave instabilities under design 
thrust oscillation conditions. 

A. Computational Modeling Tool 

The computational software used to study the tank vertical sloshing phenomenon is the commercially available 
CFD-ACE+ program, which was originally developed by CFD Research Corporation (CFDRC), and is currently 
owned and distributed by ESI [3]. CFD-ACE+ is a multi-physics and multi-disciplinary simulation tool, and is 
especially suited for vertical sloshing modeling. As thrust oscillation is a phenomenon of structure vibration induced 
by resonant burning, we will first consider the structural vibration coupled with the gas-liquid interface dynamics. 
CFD-ACE+ has a built-in coupled fluid- structure interaction module, and can be used to study the resonance of a 
tank-spring system. 

CFD-ACE+ solves the Navier-Stokes equations in a Lagrangian-Eulerian frame. The continuity and 
momentum equation can be generally written as: 

7 - lpdV + {(p(v-v g )-ds)=0 ( 1 ) 

dt V s 
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( 2 ) 


— J p(f»dV + j p(f»(v - v g )• ds = j q • ds + jS^dV 

dt V s s V 

where § are the Cartesian velocity components, v is the absolute fluid velocity, q the diffusive flux and the 
volume sources. V is the computational cell volume and S are bounding cell surfaces, v g is the grid velocity. As 
the grid is moving with time for the present fluid- structure interaction problem, a space conservation law (SCL) is 
enforced during the grid deformation, 

7 - \ dV = J v g ds (3) 

dt v 

CFD-ACE+ also features a structural dynamics module, FEMSTRESS. Its main capabilities include: 

• triangular, quad, tetrahedral, prisms, or brick elements; 

• linear or high-order isoparametric elements; 

• small or large deformation; 

• elastic or plastic stress; 

• isotropic or anisotropic materials; 

• thin to thick shell/plate elements; 

• modal analysis and eigenvalue solutions; and 

• steady and dynamic analysis. 

FEMSTRESS uses the finite element method and solves the equation of motion: 

NM+ [cKq}+ [K]{q}= M (4) 

where {q} is the displacement vector, [M] is mass matrix, [C] is the damping matrix, [K] is the stiffness matrix, 
and {F} is the force vector due to the fluid dynamic load and shear stresses. 

FEMSTRESS supports several boundary conditions specified from a Graphic User Interface (GUI): 

• specified load: pressure or concentration loads; 

• coupled pressure load from fluid; 

• specified body forces; 

• specified displacements: zero (fixed) or prescribed motion; 

• symmetry conditions; and 

• contact condition: contacting surfaces and target surfaces for elastic-rigid contact, and, elastic-elastic contact. 

Related to the present free surface, CFD-ACE+ contains a Volume of Fluid (VOF) module which is designed for 

applications involving two immiscible fluids. In the current application, the first fluid is liquid LH2, L02, or 
water, and the second fluid is gaseous H2, 02, helium or water vapor, respectively. In the VOF module, a single set 
of momentum and continuity equations is solved, but different property sets are defined for each fluid. The volume 
fraction of one phase (in this case the liquid phase) is tracked throughout the solution to determine which fluid 
occupies each computational cell at any given time. In cells containing both fluids, a special routine is used to 
locate the shape, location and normal of the interface. When there is a surface tension force, its effect is applied in 
a conservative form [4] . For time dependent simulations such as the present tank under oscillatory force, a special 
second order algorithm is used to update the volume fraction in a cell from one time step to the next. The second 
order geometric reconstruction scheme for the interface representation is employed to track the interface. This 
unique reconstruction algorithm is currently available only for structured grid. It is due to this algorithm that the 
present CFD solution is capable of capturing and maintaining a sharp interface between the phases for very long 
period of time. The time step size is determined by the local Courant-Friedrichs-Lewy (CFL) number. Some of the 
validation studies and applications to ARES I vertical and side sloshing were presented in a previous JANNAF 
paper [5]. The validation study for the tank sloshing is reported in reference [6]. 

II. Results And Discussion 

Modeling vertical sloshing under thrust oscillation is a challenging task. It involves structural dynamics 
resonance, fluid- structure coupling, and gas-liquid interface dynamics. The structural dynamics equations are 
usually written in the Lagrangian frame, whereas the fluid dynamics equations are in the Eulerian frame. Some type 
of moving mesh has to be invoked to match the moving structure and fluid grid at the structure-fluid interface. 
This study will take a building block approach: building up to the final complete 3D propellant tank model 

advancing through a series of simple, reduced dimensionality elemental verification and validation cases. First, a 
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2D axisymmetric, multi-phase model of gases, liquids, and solid structures is used to verify that CFD-ACE+ is 
capable of modeling fluid- structure interaction and system resonance. The 2D model gives fast turn around time 
and builds confidence in the accuracy of the final full 3D fluid- structure and VOF interaction study. Then, the 3D 
mode will be studied with the assumption of oscillatory body force to simulate the thrust oscillating effect. 



Figure 1. Experimental Setup and Simulation Model of the Tank and Elastic Structure 


A. 2D Self-Vibration Fluid-Structure Interaction Model 

The basic experimental set up at NASA MSFC [7] and the computational model are described in Figure 1 . The 
tank has an 8” diameter and 15” height, with a 75% initial water fill level, and is supported by elastic springs. A 
shaker acts at the bottom of the tank to provide the driving force for the “thrust oscillations”. The maximum force 
from the shaker is rather small, but it can cause thrust oscillations when resonating with the system. The natural 
frequency of the system is primarily determined by the elastic stiffness of the springs and the mass of water in the 
tank. The computational model simulating the spring-tank interaction and free surface dynamics is presented in 
Figure 1 . Here the bottom wall of the tank is modeled as an elastic solid along with the springs, and the side wall is 
treated as rigid without mass. The interface between the liquid and bottom wall is treated as a fluid-structure 
interface, where the continuity of displacement, velocity, normal stresses and shear stresses from the fluid side and 
the solid side are satisfied at all times. The fluid dynamics is modeled using Finite Volume Method (FVM) and the 
free surface is tracked by VOF. The springs are modeled as elastic structures, and are calculated by Finite Element 
Method (FEM). The computational grid and boundary conditions are shown in Figure 2. 

B. Self-Vibration Modeling 

First, the self-vibration of the system is studied by initially setting the spring at an unstressed state. At t=0, the 
gravitational force of the tank is applied to the spring, so that the self vibration starts. Since the spring is modeled 
as a solid block, the equivalent Young’s modulus has to be determined. In this study, the value of Young’s modulus 
is selected such that the system will resonate at 15Hz, the same resonant frequency as in the experiment. Figure 2 
shows the time history of the vertical displacement at a monitor point inside the solid block. During a one second 
time frame, the spring displacement shows 15 cycles, giving an oscillating rate of 15 Hz. Due to fluid dissipation, 
the amplitude of the oscillation decays, and eventually, the oscillation stops, when the spring force balances the 
mass of the tank. From an animation of Figure 2, one can clearly see the spring compression and expansion as the 
tank moves up and down. This example demonstrates that the essential physics of coupled fluid-structure 
interaction and free surface dynamics involving three phases of solid, gas and liquid are well captured. 
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Figure 2. 

Modeling of Fluid-Structure Interaction of Self-vibration of the Tank and Spring System. 


C. Forced Vibration to Simulate Thrust Oscillation. 

Next, a time dependent pressure force is applied to the bottom of the tank as shown in Figure 3. The force is 
oscillatory and at the natural frequency of the system of 15Hz. Resonance occurs, which is indicated by the 
increase of vibration amplitude of the spring plotted in Figure 3. The displacement history at the same location as in 
Figure 2 is shown in Figure 3. The amplitude increases with time, and eventually settles down as a “limit cycle”, 
where the supplied energy is balanced by viscous dissipation. The corresponding liquid-gas interface contours in 
Figure 3 show small waves that grow and eventually become standing waves. The surface dynamics is the same as 
observed in the experiments, except that the waves are in the form of concentric rings due to the 2D limitation, 
instead of droplet type of peaks and valleys in 3D. 



Shaker: P=-62.66 sin(2*PN5*t) psf 


Figure 3. Fluid-Structure Interaction for Forced Vibration. 


D. 3D Tank Modeling 

The above example clearly demonstrated the unique capability of CFD-ACE+ in capturing the fundamental 
dynamics of thrust oscillation induced surface wave, the coupled fluid- structure interaction, and multi-phases (solid, 
gas, and liquid) dynamics. In order to quantitatively compare with analytical solution and with experiment result, a 
3D computational model was generated. To save computer resources and for fast turn around time, the Thrust 
Oscillation (TO) effect was modeled as a body force, such that the tank is set in a non-inertial frame which moves 
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with the rigid tank wall. The fluid is assumed to be isothermal, and the physical properties are constant. The 
diameter of the tank is 8”, as in the experiment, and a volume grid with a ‘butterfly’ structure featuring 
approximately 60,000 cells was generated inside the tank as illustrated in Figure 4. The initial condition of the 
liquid interface is shown in Figure 4, with the liquid phase marked in red and the gas phase in blue. The upper 
boundary is open to the ambient, and the remaining boundaries are defined as walls with zero slip velocity. The TO 
effect is modeled as time dependent body force with a parameter a, which relates to the vertical body force as: 

g y = g *(l + orsin(2^3r)) (5) 

The value of a represents the g level of the oscillatory acceleration in the vertical direction. 

CFD Model 


9y 


g y = g*[l + asm(2nf*t)] 

Figure 4 3D CFD model for 8” Figure 5. Stability charts for the solution of Mathieu’s equation [2]. 
cylinder tank. 

E. Initial Linear Stability Regime. 

Natural Frequency of the Surface Wave 

The interfacial dynamics of a cylindrical tank containing a heavy liquid forced to vibrate vertically has been 
studied in 1960’s using both analytical and experimental methods. The surface dynamics is called Faraday’s wave, 
and a comprehensive review can be found in Chapter 8 of reference [2]. With the assumption of incompressible 
and irrotational flow, and with neglected gas phase effects, a neutral stability equation has been derived, which is in 
the standard form of Mathieu’s equation. The eigenmodes of Mathieu’s equation give stability charts as shown in 
Figure 5. For each eigenmode, there is a modal frequency and a modal shape function. The modal frequency is the 
natural frequency of the surface wave, and the modal shape represents the free surface wave form. 

One of the eigenmodes of the current 8” tank is shown in Figure 6(a) (Courtesy of Dr. Bernard Beard of ARES 
Corp.). The color represents the shape of the wave, with red for the peaks, and the blue for the valleys. Figure 6(a) 
is the mode of (6,4) with eigenvalue of 14.995 Hz, which is close to the driving frequency under the consideration. 
There are other modes which have natural frequency close to 15Hz for the 8” tank. 
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a) Analytical solution. (6,4) mode shape @1 4.995 HZ from (b) Standing surface wave form from CFD-ACE Simulation. 

Mathieu equation (Courtesy of Dr. Bernard Beard of Ares Combination of (4,6) mode and (3,7) mode. 

Corporation) 

Figure 6. Comparison of surface wave form and wavelength between analytical solution and CFD 
simulations. 

Surface Wave Response to the Excitation Forcing Frequency 

A simpler version of Figure 5 is shown in Figure 7(a) [2]. Here, cOn represents one of the natural frequencies 
close to the forcing frequency. The horizontal axis represents the driving frequency parameter. The vertical axis 
in Figure 7(a) is the wave amplitude with the curves presenting the neutral stability. For any state above the curve, 
the surface wave is unstable, and for the state below the curve, the wave is stable. When the excitation frequency, 
which is the tank-spring resonant frequency for the current system, coincides with one of modal frequencies 
(frequency parameter is equal to 1.0 in Figure 7(a)), a small amplitude forcing will lead to the surface wave growth. 
However, when the excitation frequency does not coincide with the natural frequency of the surface wave, the free 
surface will become unstable only when the amplitude of the excitation is larger than a certain value. 

Onset Frequency and Forcing Amplitude Validation 

To study this amplitude dependent instability, numerical simulations were carried out using various a values 
from 0.1 up to 2.0 at 15Hz (see Equation (5) for the definition of a). All simulations started with the quiescent no 
motion state. The interface wave amplitudes were recorded from the CFD simulations and are shown in Figure 
7(b). At smaller values of the excitation force, when the a value is below 0.8, there is no sloshing motion at all. 
The surface of the liquid-gas is essentially flat, and the liquid column acts as a rigid body and vibrates up and down. 
The pressure field stratifies instantaneously with the time dependent body force variation. However, when the 
oscillatory acceleration g level (a) is above 0.8, surface instability sets in, and a pattern of standing waves can be 
observed on the free surface. The wave amplitude is found to increase with increasing a as seen in Figure 7(b). 
This evolution of dynamics is just as predicted from the nonlinear solution of Mathieu’ s equation under large 
amplitude shown in Figure 7(a). If the driving frequency is off the natural frequency of the surface wave, only 
when the amplitude of the driving force is above the neutral value, will the instability set in. 

Wavelength and Wave Form Validation 

To compare the surface wave shape, the CFD solution is given in Figure 6(b), along with the analytical solution 
of the modal shape expressed in terms of Bessel functions. It should be noted that there are several natural modes 
for this 8” diameter tank with frequency close to 15Hz. Shown in Figure 6(a) is just one of these frequencies, with 
the mode (6, 4). The CFD simulation shows very similar wavelength and wave pattern, but it may well be a 
combination of the (4, 6) and (3, 7) modes. One should also notice that the waveform of Figure 6(b) is in large 
amplitude regime. Regardless, agreement between the analytical solution and CFD prediction is fairly good. To 
make further comparison with experiment, the photographs from the test along with surface shape from CFD 
simulation are given in Figure 8. Again, the agreement is good in terms of wavelength and wave form. 
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(a) Analytical Solution from Mathieu’s Equation [2] (b) Present CFD Solution Analytical Solution 

Figure 7. Stability curve from CFD solution and from analytical solution. 


Surface Wave Response Frequency Validation 

One of the interesting features of vertical sloshing is that the wave frequency, or the response frequency of 
surface wave, usually occurs at exactly one-half that of the tank motion. Although in some cases the liquid 
frequency is equal to the forcing frequency and in other cases is greater than the forcing frequency. This is very 
different from side sloshing where the frequency of the free surface to the transverse excitation corresponds 
precisely to the forcing frequency. This apparently anomalous behavior for vertical excitation was first noticed by 
Faraday [2] during a series of investigations of vibrating plates covered by a thin layer of water. The liquid 
frequency was recorded in the present CFD simulation, and it was found that the liquid oscillation was indeed at the 
one-half subharmonic of the forcing motion. 


Validation Summary for the Linear Stability Regime 

These simulation results again show the capability of CFD in capturing the fundamental physics in terms of a) 
onset frequency; b) surface wave form; c) driving forcing amplitude; and d) surface wave response frequency. 



Excellent agreement in wavelength and wave form 


Experimentation at MSFC 

(Courtesy of Dr. John Townsend, EV31) Present Modeling 


Figure 8. Comparison of surface wave form and wavelength between experiment testing and CFD 

simulations 
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F. Non-Linear Regime 

The above test case shows that CFD-ACE+ is capable of determining the transition characteristics and post- 
transition behavior of Faraday’s waves (not possible by linear analytical solutions) for smaller amplitude of the 
forcing frequency, i.e. transition from solid body motion to standing wave motion. The experimental work of 
Goodridge et al. [8] took the above Faraday’s waves further into the nonlinear regime, where the unbroken surface 
forms and ejects droplets under parametric excitation. They found that with increasing forcing, an initially flat 
surface wave state will evolve into a periodic standing wave state. This then will change into an aperiodic state, 
which, at sufficiently large forcing, will produce large-amplitude waves in the form of jets (or spikes) which will 
break into droplets. The purpose of this study is to investigate if a CFD tool can model the observed droplet 
ejection from the surface wave due to TO at a high amplitude of the forcing, and to determine the grid resolution 
requirement to resolve the droplets from the first principle. 



g y = g*[l + asin(2;^*0] 


Experimental 3D Cylinder Model of 
Goodridge et ai. [Phys. Rev., vol. 56, pp. 
472-476, 1997] 



40 cells 
0.087cell 


_40 cells, 0.02'7cell 
40 cells, 0.02"/cell 


40 cells 
0.087cell 


CFD 2D Cartesian Model 


Figure 9. Experiment and simulation model to study the droplet ejection of Faraday’s waves under large 

driving force. 


The experimental configuration of Goodridge et al. [8] is shown in Figure 9, where a cylindrical plastic 
container of 20 cm (7.677”) diameter with initial liquid filled to 10 cm (or 3.93“) is mounted on the armature of a 
electromagnetic shaker. The liquids used in the experiment were water, water-glycerin, and ethanol, representing 
different viscosities. Our study was first done using a 2D model, where 400 grid cells were spaced evenly along the 
horizontal direction, and variable grid sizes were used along the vertical direction with a denser mesh near the 
interface to ensure that sufficient grid resolution was applied to capture the high frequency small wavelength waves. 

The experimental work of the Goodridge et al. [8] covered the forcing frequency from 20Hz up to 80Hz. Under 
high frequency forcing, the response wavelength of the natural mode decreases with the increase in the forcing 
frequency. Depending on the wavelength, the restoring energy from ejection could be either surface tension or 
gravitational force. At smaller wavelengths, surface tension can be dominant. To determine the critical 
wavelength between the two forces, one can derive a dispersion relationship as shown in [10] for the two- 
dimensional gravitational wave of Figure 10 as: 

2 , &k 3 

co=gk + (6) 

p 


where co is the angular frequency, g is the gravitational acceleration, k is the wave number, a is the surface 
tension, and p is the liquid density. With k=2n/X, X as the wavelength, one can find the frequency as: 


/ = 


8 , 

2nX pX 3 


( 7 ) 
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Figure 10. Surface gravitational wave on liquid and gas interface 


It can be seen that the wave frequency depends on surface tension, gravity acceleration, and wavelength, 
critical wavelength can be found when the surface tension contribution is equal to that of gravity: 

g _ 2 no 

2jt1 o'/ L 3 


( 8 ) 


A„ = 2n „ I — (9) 

PS 

At ground conditions (sea-level gravity) with water properties, one can find that /l = 0.017m . Thus, when 


the wavelength of the free surface is less than 0.017m, or 0.67”, the surface tension force becomes the main 
restoring force, and capillary waves will disintegrate. Based on the critical wavelength, we can also find a critical 
frequency of 13.8Hz. When the forcing frequency is below 13.8Hz, the restoring force is the gravitational force, 
whereas when the forcing frequency is larger than 13.8Hz, the surface tension force becomes important. In the 
following discussion, In the following discussion reference will be made to a critical wavelength greater than 0.67” 
and corresponding critical frequency of 13.8 Hz. 


Driving Frequency at 20Hz 

First we consider the case when the driving frequency is at 20Hz. At this frequency, the resulting surface wave 
dynamics are dominated by the surface tension. As the resulting wavelength is below 1.7cm, the container 
boundaries (which are 7.67 cm apart) should have less influence on the results. Representative results are given in 
Figure 11. The experimentally determined threshold value of acceleration from reference [8] is shown at the lower 
right corner of the figure. For the case with water as the liquid, the ejection threshold of alpha is slightly larger than 
1 .0 from the test data. 

Two simulations are shown in Figure 11, one using a of 1.0 and one of 2.0. For the case with oc=1.0, the 
standing wave is visible and there are no droplets stripped from the surface. The surface integrity is intact. The 
standing wave leads to a surface area increase which may potentially increase heat transfer between liquid and gas 
phase, but this increase is not significant. However, when the forcing amplitude increases to oc=2.0, as shown in the 
same figure, the surface standing waves disintegrate and form a large amount of droplets from the surface. In the 
animation (not shown in this paper, but in the presentation), the droplets form a dense spray, and they remerge with 
the liquid surface. Due to high initial break-up velocity, the droplets could enter into the region beyond the 
interface neighborhood. For propellant tanks, these droplets can greatly enhance heat transfer between liquid and 
gas phases. Not only will they carry heat away from the gas phase near the free surface due to phase change, but also 
they carry heat away from the warmer gas well within the ullage region. This type of droplet ejection can lead to 
potential ullage collapse. 

An experimental photograph from Reference [8] is used for comparison with the CFD simulation in Figure 11. 
Again, good agreement is evident in terms of wavelength, wave pattern, and spawned droplets. Favorable 
comparison is also found in the threshold acceleration value shown in Figure 11, where oc=l is below the threshold 
and oc=2 is above the threshold. As the droplet sizes are all smaller than 0.67”, the surface disintegration is surface 
tension dominated. It should be pointed out that the smallest droplet size that can be captured in the simulation is 
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about the grid cell size of 0.02”. Objects smaller than the grid cell cannot be resolved numerically with the present 
phase interface tracking algorithm. 



f=20Hz, a =2.0 

Droplet ejection from experiment and simulation 


□ 85% Glycerin- Water A 

A 84% Glycerin- Water D 
100 [_+ <0% Glycerin- Water g A 
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f=20Hz, a =1.0 
No droplet ejection from experiment and simulation 


Good agreement between present CFD model and experiment. 
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Threshold Acceleration from Goodridge et at. (1997) 
Figure 11. Droplet ejection modeling and comparison with experiment at 20Hz using water 


G. Driving Frequency at 60Hz 

When the forcing frequency increases, the disintegration process is suppressed at the same forcing amplitude of 
a=2. This is evident in Figure 12. With a=2, there is no droplet ejection from the surface, rather a standing wave 
is observed. The CFD simulation is consistent with the experimental observation, where the threshold increases with 
increasing forcing frequency (see the experimental curve in the lower right hand side of Figure 12). One can also see 
that the wavelengths have reduced in comparison with those at the case of f=20Hz. However, for a large driving 
force, a large amount of small droplets can be seen stripped from the surface. Comparison of droplet sizes and 
surface wavelength between CFD and experiment is again fairly good. The threshold value is bounded from the 
present CFD simulations and its value is between 2 < a < 4. 
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Droplets being ejected 
from surface. The 
smallest size: 0.02" 
diameter (the grid cell 
size) 


f=60Hz, a =4.0 

Droplet ejections from experiment and simulation. Shorter surface wavelength, 
smaller droplet size compared to those at f=20Hz. 
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Good agreement between present CFD model and experiment. 
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Threshold Acceleration from Goodridge et al.(1997) 


Figure 12. Droplet ejection modeling and comparison with experiment at 60Hz using water 


H. Viscous Effects in Droplet Ejection Waves 

The above comparisons were made for air and water as the fluids. The experiments of Goodridge et al. [8] 
using different fluids found that low-viscosity fluids have threshold accelerations which depend on only surface 
tension and forcing frequency, while high viscosity fluids have thresholds which depend on only viscosity and 
forcing frequency. This effect is shown in Figure 11 with various fluids. In order to simulate the experimental 
observation, CFD simulations were made using a glycerin- water mixture with a viscosity 43 times higher than 
water. The simulated instantaneous surface wave and the corresponding experimental photograph are shown in 
Figure 13. In contrast to water where the wave breaks up into droplets immediately, here several long filament-like 
peaks are visible before break off occurs. The high liquid viscosity enables the peaks to achieve their filament 
structure. The comparison with experiment in terms of the wave structure is rather remarkable as seen from Figure 
13. The higher viscosity can stabilize the droplet ejection, as shown in Figure 14, where the same driving 
frequency and acceleration level oc=2 are shown for water and glycerin-water. The small droplets observed in water 
are suppressed in glycerin- water. 
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Effect of Viscosity: Water with lower viscosity 
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f=20Hz, a =2.0 


Effect of Viscosity: Glycerin -Water with higher viscosity 



f=iOHz, a=3.0,p=43xp^ ter 

Different from water, where wave breaks up into droplets immediately. Here 
several long filament like peaks are visible before break off occurs. 

Figure 13. Droplet ejection modeling and comparison with experiment at 60Hz using Glycerin-water 


g> = g*[l + asin(2^*0] f=20Hz 



Water, a=2.0 


Glycerin-Water, cl=2.0 


Glycerin-Water, a=3.0 


Figure 14. Characteristics of droplet formation at different acceleration level and for different fluids. 


III. Summary And Conclusions 

The multi-phase fluid interaction simulation capability of the CFD software program, CFD-ACE+, has been 
demonstrated and validated on several key phenomena of liquid surface instability subject to vertical thrust 
oscillation. The present study showed that CFD modeling can accurately capture all the transition physics from solid 
body motion, to standing surface wave, to droplet ejection at the liquid-gas interface. The present simulation also 
demonstrated trends consistent with experiments in the variation of physical properties, geometry, and input forcing 
function. It was found that thrust oscillations cause increased phase exchange area at the interface due to droplet 
formation, and enhanced heat transfer rate between the liquid and gas phases as the ejected droplets travel well into 
the warmer gas region. The first-principle based CFD simulation can capture the detailed droplet ejection as long the 
grid size is not larger than the droplet diameter itself. 

The present effort also builds a solid foundation for the practical applications of CFD technology to address the 
space exploration systems involving control and management of cryogenic propellant tanks using thermodynamics 
vent, and to the advanced propulsion systems involving liquid jet breaks-up and atomization. 
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